set.seed(2423045)
rm(list=ls())
library(Matching)

source("AutoCluster3.R")

load("rep.sourcecode.nore74psid1.RData")

X <-   as.matrix(cbind( foo$age,
                        foo$educ,
                        foo$black,
                        foo$hispan,
                        foo$married,
                        foo$nodegree,
                        I(foo$re75/1000)))# ,
                        #I(as.real(foo$re75 == 0))))
                 
BalanceMat <- as.matrix(cbind(
                     foo$age,
                     foo$educ,
                     foo$black,
                     foo$hispan,
                     foo$married,
                     foo$nodegree,
                     I(foo$re75/1000),
                     I((foo$re75/1000)^2),
                         
                     I((foo$age)^2),
                     I((foo$educ)^2),
                     I((foo$black)^2),
                     I((foo$hispan)^2),
                     I((foo$married)^2),
                     I((foo$nodegree)^2),
                         
                     I(foo$age*foo$educ),
                     I(foo$age*foo$black),
                     I(foo$age*foo$hispan),
                     I(foo$age*foo$married),    
                     I(foo$age*foo$nodegree),
                     I(foo$age*(foo$re75/1000)),

                     I(foo$educ*foo$black),
                     I(foo$educ*foo$hispan),
                     I(foo$educ*foo$married),    
                     I(foo$educ*foo$nodegree),
                     I(foo$educ*(foo$re75/1000)),

#                     I(foo$black*foo$hispan), This is always zero!
                     I(foo$black*foo$married),
                     I(foo$black*foo$nodegree),
                     I(foo$black*(foo$re75/1000)),

                     I(foo$hispan*foo$married),
                     I(foo$hispan*foo$nodegree),
                     I(foo$hispan*(foo$re75/1000)),

                     I(foo$married*foo$nodegree),
                     I(foo$married*(foo$re75/1000)),
                                                
                     I(foo$nodegree*(foo$re75/1000)) 
                         ) )


# From Review of Econ and Statistics
model.A = foo$treat~I(foo$age) + I(foo$age^2) + I(foo$age^3) + I(foo$educ) + I(foo$educ^2) + I(foo$married) + I(foo$nodegree) + I(foo$black) + I(foo$hispan) + I(foo$re75) + I(foo$re75 == 0)

pscores.A <- glm(model.A, family = binomial)

X2 <- cbind(X,BalanceMat[,8:ncol(BalanceMat)]) #first 7 spots of BalanceMat are already in X
for (i in 1:ncol(X2))
  {
    lm1 = lm(X2[,i]~pscores.A$linear.pred)
    X2[,i] = lm1$residual
   }

# VARIABLES WE MATCH ON BEGIN WITH PROPENSITY SCORE
orthoX2.plus.pscore = cbind(pscores.A$linear.pred, X2)
orthoX2.plus.pscore[,1] <- orthoX2.plus.pscore[,1] -  mean(orthoX2.plus.pscore[,1])


# NORMALIZE ALL X COVARS BY STANDARD DEVIATION
#for (i in 1:(ncol(X)+1))
for (i in 1:ncol(X2))
  {
    orthoX2.plus.pscore[,i] <-
      orthoX2.plus.pscore[,i]/sqrt(var(orthoX2.plus.pscore[,i]))
  }


#sv from psid1.nore74.gm1.triton1.Rout
sv <- c(9.083644e+02, 8.042862e+00, 4.405689e-01, 3.514126e+00, 2.135189e+00,
        8.834534e+02, 5.168435e+02, 4.553718e+02)
sv <- c(sv,rep(0,ncol(orthoX2.plus.pscore)-length(sv)))

#orthoX2.plus.pscore <- orthoX2.plus.pscore[,1:(8+cnew)]

for (s in 1:100)
  {
    cat("*S:",s,"**************************************************************\n")
    cl <- NCPUS()
    GM.out <-  GenMatch(Tr = foo$treat,
                        X = orthoX2.plus.pscore,
                        starting.values=sv,
                        BalanceMatrix = BalanceMat,
                        pop.size = 10000,
                        max.generations = 500,
	                hard.generation.limit=TRUE,
                        wait.generations = 50,
                        cluster=cl)
    stopCluster(cl)

    sv=diag(GM.out$Weight.matrix)
    cat("Starting Values:\n")
    write.csv(sv,eol=", ", row.names=FALSE)

    mout <- Match(Y=foo$re78,
                  Tr = foo$treat,
                  X = orthoX2.plus.pscore,
                  Weight.matrix=GM.out)
    summary(mout)
  }


